
% code for averaging of SPMs wrt coherence
% Reference : Stodt J.A, 1983.
% Uses hard limiter (1 &0 )
% NOT tested

% Latest date 25.03.2003


function[SPM,ProcDef] = cohspm(SPMatrix,ProcDef,m),
if ProcDef.overlap ==1,  % Overlap of windows cause reduction in dof for FC's
	v=0.818;					 % A correction to this effect is done here
else,
 	v=1;
end;

A = size(SPMatrix);
trsh = ProcDef.CohTr;

if trsh < 0 | trsh > 1,
   fprintf('cohspm-> Warning trsh is reset to 1.0\n');
   fprintf('         No change from this module\n'); 
   trsh=1;
end;

%CHxEy = ocoh(SPMatrix,5,1); % returns ordinary coherency 
%CHyEx = ocoh(SPMatrix,4,2); % between H & E fields; 2D array [nfrq len]
CHxEy = pcoh(SPMatrix,'Ex'); % returns ordinary coherency 
%CHyEx = pcoh(SPMatrix,'Ey'); % between H & E fields; 2D array [nfrq len]
%C = sqrt(CHxEy.*CHyEx);
C = CHxEy;


if m == 1, % Preferential stacking based on coherence
  	for i = 1:A(1),
		h = sort(C(i,:));
		tre = h(floor((0.8)*A(2)));
		k = find(C(i,:)<tre);
		ProcDef.Weights(i,k) = 0;
	end;
elseif m == 2, % Coherence threshold stacking
   for i = 1:A(1),
      k = find(C(i,:)<trsh);
      if A(2)-length(k) > 4,
         ProcDef.Weights(i,k) = 0;
      end;
   end;
end;

for i = 1:A(1),
   SPMat(:,:,:) = SPMatrix(i,:,:,:); % all the events for one frequency
   SPM.spectra(i).data(:,:) = WeightSPM(SPMat,ProcDef,i);
   SPM.spectra(i).freq = ProcDef.TLFreq1(i);
   SPM.spectra(i).dof.avgt = floor(sum(ProcDef.Weights(i,:))*v); 
   SPM.spectra(i).dof.avgf = ProcDef.TLRad1(i)*2; % AVGF basically results from n of parzen harmonics

end;
 SPM.head.site = ProcDef.file1;
